PASJ: Publ. Astron. Soc. Japan , -?? (2008) 



Galaxy Formation from a Low-Spin Density Perturbation 

in a CDM Universe 



Or 



> 
(N 
(N 

o 

On 
On 

S3 

Oh' 

6 

CO 



X 



Daisuke Kawata 
Astronomical Institute, Tohoku University, Sendai, Miyagi 980-8578 
E-mail( DK): kawata@astr. tohoku. ac.jp 

(Received ; accepted ) 
Abstract 

In order to understand the formation process of elliptical galaxies which are not rotationally supported, 
we have carried out numerical simulations of the galaxy formation from the density perturbation with a 
rotation corresponding to a small spin parameter. The three-dimensional TREE N-Body/SPH simulation 
code used in this paper includes the dark matter and gas dynamics, radiative cooling, star formation, super- 
nova feedback, and metal enrichment. The initial condition is a slowly rotating, top-hat over-dense sphere 
on which the perturbations expected in a cold dark matter (CDM) universe are superposed. By means of 
the stellar population synthesis, we calculated the surface brightness profile, the metallicity distribution, 
and the photometric properties of the end-product, and found that these properties quantitatively agree 
with the observed properties of bright elliptical galaxies. Thus, we conclude that, in a CDM universe, the 
proto-galaxy which has a spin-parameter as small as 0.02 evolves into an elliptical galaxy. 
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1. Introduction 

The new generation of large ground-based telescopes 
and the Hubble Space Telescope allow us to observe 
structures and properties of galaxies in unprecedented 
detail. In order to clarify the history of galaxy formation 
and evolution based on the observed features of galax- 
ies, the numerical simulations are required which can be 
compared with these observations directly. The numer- 
ical models of galaxy formation have progressed greatly 
since it was first developed by a series of Larson's pa- 
pers. Larson (1969) modeled the gas dissipation and the 
chemical evolution in galaxy formation and reproduced 
main properties of the observed elliptical (Larson 1974a, 
1974b, 1975) and disk (Larson 1976) galaxies success- 
fully. However, he did not take into account the details of 
the hydrodynamics and dissipation of the gas, and only 
considered the axisymmetric system. Carlberg (1984a, 
1984b) and Carlberg et al. (1986) simulated the dissi- 
pational formation of spheroidal galaxies using N-body 
model. They succeeded in reproducing elliptical galaxies 
with a de Vaucouleurs surface density profile. Although 
they improved Larson's approach, some problems still 
remained. For instance, the treatment of the gas compo- 
nent and of the star formation process was phenomeno- 
logical, because they introduced various parameters, such 
as the collision rate of the clouds, the rates of dissipation 
and star formation after encounter between gas particles, 
and the strength of the kinetic feedback to the surround- 



ing gas caused by supernovae of new born stars. Further- 
more, their simulations were performed in dimensionless 
units, so that there was no physically convincing way of 
scaling their results of simulations to dimensional units 
and comparing them with the observational data. 

Katz & Gunn (1991) presented the modeling galaxy 
formation using the general-purpose code for evolv- 
ing self-gravitating fluids in three dimensions, called 
TreeSPH (Hernquist, Katz 1989; Katz et al. 1996). 
TreeSPH is a fully Lagrangian code to combine smoothed 
particle hydrodynamics (SPH; Lucy 1977; Gingold, Mon- 
aghan 1977) with a hierarchical tree algorithm for compu- 
tations of gravitational forces (Barnes, Hut 1986). Since 
in this method the thermodynamic state of the gas can 
be calculated, the dissipational effects of the gas are in- 
cluded as the radiative cooling following the standard 
cooling curves. Thus, simulations based on this numeri- 
cal model can lead to physical values which can be com- 
pared with the observational data quantitatively Fur- 
thermore, the TreeSPH simulation was presented which 
includes the process of the star formation (Katz 1992) 
and of metal enrichment due to supernova (Steinmetz, 
Muller 1994, 1995). 

Using these numerical model, they (Katz, Gunn 1992; 
Katz 1992; Steinmetz, Muller 1994, 1995) investigated 
the collapse of an isolated constant-density perturbation 
which initially follows a Hubblc-flow expansion and has 
a solid-body rotation, since external tidal fields are not 
included. The top-hat density perturbation consisted 
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of dark and baryonic matter, and included small scale 
perturbations assuming power spectrum index —2.5. In 
their model, the amount of initial solid-body rotation was 
specified by the dimensionless spin parameter, A (Peebles 
1971). They studied only the case with A = 0.08, and 
seemed to succeed in forming a system which had some 
similar properties to an observed disk galaxy. However, 
some groups (e.g., Barnes, Efstathiou 1987; Warren et 
al. 1992) concluded that the spin parameter of a virial- 
ized dark matter in a cold dark matter (CDM) universe 
falls on a range of 0.02 — 0.11 with the median value of 
0.05. The recent analytical studies also derived the dis- 
tribution of spin parameter with the similar spread and 
median value (e.g., Steinmetz, Bartelmann 1995; Eisen- 
stein, Loeb 1995), in which A = 0.08 is near a high end 
value. By this reason, it is worth while investigating the 
properties of the system formed from the initial density 
perturbation with a small value of A. 

Therefore, we calculate the dynamical and chemical 
evolution of density perturbations with the small spin 
parameter, A = 0.02. The present simulation follows a 
semi-cosmological manner (Katz, Gunn 1991) in order 
to specify initial density and velocity fields, and adopts 
the method of numerical simulation used by Katz (1992) 
and Steinmetz & Miiller (1994, 1995). Using stellar pop- 
ulation synthesis (Kodama, Arimoto 1997), we can di- 
rectly compare the results of simulations with the ob- 
served properties of galaxies, such as the surface bright- 
ness profile, the half-light radius, the metallicity gradient, 
and the color gradient. 

According to Heavens & Peacock's (1988) analytical 
estimation, high density peak tends to have a small value 
of A. This reason is that higher peaks have a shorter col- 
lapse time; thus they have less time to get spun up by 
the environmental tidal force. On the other hand, bright 
elliptical galaxy is considered as the system which col- 
lapsed at high redshift and has a small rotation. It is a 
natural supposition that ellipticals were formed from a 
high density peak (e.g., Blumenthal et al. 1984). Hence, 
we focus on comparison the results of numerical simula- 
tion with the observed properties of elliptical galaxies. 

Section 2 presents the method of numerical simula- 
tions, and the details of the top-hat model of galaxy for- 
mation are described in section 3. Section 4 presents the 
results of numerical simulations. Summary and discus- 
sion are described in section 5. 

2. The Code 

In our numerical simulation, the dynamics of the dark 
mater and stars is calculated by N-body scheme, and the 
gas component is modeled using the SPH. The code also 
includes the processes of the gas cooling and of the star 
formation. Many authors (e.g., Hernquist, Katz 1989; 
Navarro, White 1993; Katz et al. 1996; Steinmetz 1996; 



Carraro et al. 1998) have already described the details 
of numerical model for galaxy formation which employs 
the SPH method. Our code is based on these previous 
modelings essentially and is the revised version of the 
SPH code in Kawata, Hanami (1998) for the application 
to galaxy formation. We describe details of our SPH code 
and the model of the cooling and star formation in the 
following. 

2.1. The SPH Code 

In the SPH method, the densities and other local quan- 
tities in the gas component are computed by a kernel 
estimation. We employ a spherically symmetric spline 
kernel (Steinmetz 1996), 

W(r,h) = 

l-6{r/h) 2 +6{r/h) 3 if < r/h < 1/2, m 



X 



2[1 - (r/h)Y 




if 1/2 < r/h < 1, 
otherwise. 



Each particle has its own smoothing length, hi , chosen so 
that a fixed number of neighbor particles are contained 
within its smoothing length (e.g., Hernquist, Katz 1989). 
Using this kernel, the smoothed density associated with 
the i-th particle is given by 



where = \x 



(2) 

hj)/2 is the pair- 



Xj | and hij = (hi 
averaged smoothing length. 

The acceleration of the i-th particle is determined by 
Euler's equation, 



at pi 



Prise) - Vi$, 



(3) 



where <!>, P, and P v isc are the gravitational potential, the 
pressure, and the artificial viscous pressure, respectively. 
The pressure-gradient term can be written as (e.g., Stein- 
metz 1996) 



±Vi(P + P v 



= Ej (f + J + Qii) ViWixij, h^), 



where 
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if Xij ■ Vij < 0, 
otherwise, 
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Here, we define Xj 



Xi — x+ and Vt 



(4) 



(5) 



(6) 



In 



this paper, we chose parameters (a, (3) = (0.5,1.0) in 
equation (@) and j] = 0.05/iy in equation (||) so that 
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the features of the one-dimensional shock tube can be 
well reproduced. We use the corrected artificial viscosity, 
which was proposed by Navarro & Steinmetz (1997): 

Si + n 



Qij — 

fi = 



Qi 



l<v 



(7) 



|(V • v) t \ + |(V x v)i\ + 0.0002« M /V 

where v s is the sound velocity. The velocity divergence, 
(V • v)i, and rotation, (V x v) i x , of the i-th particle are 
calculated by 

1 



(V 



— y ^rrijVij ■ S/iW(xij,hij), 
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i^i, z W{x i:j , h i3 )} 



(8) 



(9) 



where Vij tX = Vi lX — Vj )X . The gravitational-force term of 
the i-th particle is the summation of the contributions 
from dark matter, gas, and star particles, 



= - G E 



(r 2 - 



.)3/2 



(10) 



where Sij = (si + Sj)/2 is the mean softening length of i 
and j-th particles. The gravitational force for each par- 
ticle is computed using the tree method (e.g., Barnes, 
Hut 1986; Pealzner, Gibbon 1996) with expansions to the 
quadrupole order and the tolerance parameter — 0.8. 

The evolution of the internal energy, itj, of the i-th 
particle is determined by the energy equation (e.g., Stein- 
metz 1996), 



dt ~ pi 



+ 2 



in i 



V,U>' (/ .//, ; i f 



7~ti 



(11) 



where A/p and 7i are the cooling rate and the heating 
rate per unit mass, respectively. In this paper, we con- 
sider the cooling which is due to the radiative process 
of the gas having a specified metallicity and the heating 
which is caused by the feedback of the massive stars (see 
following sections for detail). The pressure, P, of each 
particle can be obtained by the equation of state for an 
ideal gas, 



P = ( 7 - l)pu, 



(12) 



where 7 = 5/3 for a mono-atomic gas. 

The evolution of the smoothing length, h, is deter- 
mined so that each particle has a roughly constant num- 
ber of neighbors, i.e., the particles located within h from 
that particle. Thus, the smoothing length, h n+1 , at a 
time step n + 1 is determined by the smoothing length, 
h n , and the number of neighbors, N n , at the previous 
time step n; 

1/3" 

I + ' • - I , (13) 



h n+l =h n_ 
2 



where N s is a parameter (Hernquist, Katz, 1989). Our 
one-dimensional shock-tube tests have revealed that the 
acceptable range of N s is [34, 68] in the three dimensional 
space. In this paper, we adopt N s = 40. 

In order to update the velocity and the position of each 
particle, equation (||) is integrated using a leap-frog algo- 
rithm with individual time steps (Hernquist, Katz 1989; 
Makino 1991). The time step of the i-th particle is cho- 
sen to be Ati = min(0.3AtcFL,i, O.lAtf^), where AtcFL 
is determined by the Courant-Friedrichs-Lewy condition 
as 

AtcFL.i = — y, ; rr, (14) 

v St i + 1.2(av S: i + pmaxjl^ij]) 

and Aif is determined by the requirement that the force 
should not change too much from one time step to the 
next time step, leading to the specification 



dt 



(15) 



In a dense region, the smoothing length may become 
shorter than the softening length, requiring a very small 
time step. In order to avoid this, we set the lower 
limit of the smoothing length to be h mnl = e/2. For 
the collisionlcss particles which represent dark mat- 
ter and stars, the time step is determined by At = 
mm[0.16{s/v),0A{s/\dv/dt\) 1 / 2 }. 

These time steps do not always satisfy the constraints 
from the time scale of the internal energy evolution, be- 
cause the time scale of the radiative cooling can be very 
small. Thus, when the evolution of the internal energy 
of each particle is integrated, we do not use the individ- 
ual time step scheme but the integration is carried out 
every time step which is the minimum among all the par- 
ticles. Moreover, the energy equation is integrated using 
a semi-implicit method (Hernquist, Katz 1989). In the 
implementation, the radiative cooling is damped by the 
following equation in order to ensure numerical stability 
(Katz, Gunn 1991), 



t \ damped 



a(du/dt) iad 



y/a 2 + [(du/dt) 

u ( du\ 
0.5— + — . 

A* Wad 



(16) 



Here, (du/dt)^ is the change in the thermal energy due 
to the adiabatic compression or expansion of the gas ex- 
cluding the contribution from the undamped radiative 
cooling, (du/dt) Iad . 

Our code has been vectorized and parallelized using 
the method of Makino (1990). The parallelization is 
done only in calculating gravitational forces and search- 
ing neighbors. Since these calculations occupy most of 
the computational time, this parallelization is very ben- 
eficial. The MPI library allows our code to run on any 
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computers from the serial to the parallel one, and we can 
obtain the same results on different computers. 

2.2. The Radiative Cooling 

We employ the cooling model in Theis et al. (f992) 
who approximated the cooling rate calculated by Dal- 
garno & McCray (f972) and Bohringer & Hensler(f 989) 
as follows; 



A(T, Z, n H ) = A (Z)T ro ^n| erg cm" 3 s -1 , 



(17) 



where Z and riH are the metallicity and the hydrogen 
number density of the gas, respectively. A (Z) and m(Z) 
are defined for different temperature ranges as follows; 
at f 4 K < T < 10 5 K, 

m(Z) = 30Z, 
A (Z) = 4x i -22.i-4m(z ); 



at T > fO 5 K, 
m(Z) = 



2.b + 7y/Z 



5 - log(1.48 x lO 11 ^ 11 + 10 6 ) ' 

For temperatures above 5 x 10 4 K, bremsstrahlung be- 
comes important, which is expressed here as: 

A br = 1.4 x 10- 27 (/ ion ) 2 n o n ion T°- 5 ( 5ff ), (18) 

where (g^) « 1.35 is an averaged Gaunt factor and, 

f — ^c/^ion i 

n- - n » (l 3 Y X lz 
Ulon ~ 1-Y-Z I, 4 18 Z 

n c = nH U-ly-lz 
1-Y-Z \ 2 9 

Here, Y denotes the mass fraction of helium and we set 
Y = 0.24. We do not consider the cooling below fO 4 K, 
because our simulation does not have enough mass resolu- 
tion for cooling processes in such low temperatures (e.g., 
H2 cooling) and because the cooling rate declines rapidly 
below f 4 K. Furthermore, the heat sources such as the 
UV background radiation, which are ignored in this pa- 
per, could warm the gas components to this temperature 
easily. Thus, we set the lower limit of the temperature 
to be T lim = fO 4 K (Katz, Gunn f99f). 

2.3. The Star Formation 

We model the star formation using a method similar to 
that of Katz (1992) and Katz et al. (f996). We assume 
that the star formation occurs where the gas density is 
larger than a critical density, p cr i t = 7 x fCT 26 g cm~ 3 , 
and where the gas velocity field is convergent, V • Vi < 



0. When a gas particle is eligible to form stars, its star 
formation rate (SFR) is, 

dp* _ dp g 
~dt ~ 



dt 



t, 



(19) 



where c» is a dimcnsionless SFR parameter and t g = 
^3ir/16Gp g is the dynamical time, which is longer than 
the cooling time scale in the region eligible to form stars. 
This formula corresponds to the Schmidt law that SFR 
is proportional to p g 5 . In this paper, we set c* = 1. 



Equation (|19|) implies that the probability, , which one 
gas particle forms stars during a discrete time step, At, 



P* 



I — exp(— c^At/tg). 



(20) 



In order to avoid making an untolerably large number 
of new star particles with different masses, we set the 
whole gas particle transforming to one star particle with 
this probability. After the particle changed to the star, 
it behaves as a collisionless particle dynamically. 

We take into account the energy feedback to the sur- 
rounding gas by supernovae. We assume that each super- 
nova yields thermal energy of 10 51 ergs. For simplifica- 
tion, we assume that each massive star (> 8Af Q ) explodes 
within the simulation time step in which it was born (in- 
stantaneous recycling). In our simulation, the feedback 
also adds the ejected mass to the surrounding gas. In 
this process, a portion of metals produced in the stars is 
also returned to the gas, leading to chemical enrichment. 
We assume that an m M star yields (0.357to— 2.2) Mq 
of metal (Maeder 1987; Steinmetz, Muller 1994). The 
ejected mass is (m — 1.4) Mq, where the remnant mass 
is assumed to be 1.4 M©. The feedback energy, ^sn, 
the ejected mass, Msn, and the ejected metals, AfzgN, 
can be obtained from the Salpeter IMF (x — 1.35) with 
the lower mass limit of 0.1 M Q and the upper mass 
limit of 60 M©. This IMF leads to M SN ~ 0.122 M & , 
E SN - 0.007310 x 10 51 ergs, and Afz SN - 0.02754 Mq, 
for each 1 Mq star born. These mass, energy, and metals 
are smoothed over the neighbor particles using the SPH 
smoothing algorithm. For example, when the i-th parti- 
cle changes from the gas to a star, the increment of the 
mass of the j-th neighbor particle due to explosion of the 
new star is given by 



AM 



SN,j 



— —MsN,i, 
Pg,* 



where 



(21) 



(22) 



p g ,i = {Pg(xi)) =^2m j W(r ij ,h ij ). 
The Model 



Using the above code, we calculate a semi-cosmological 
model as follows. We consider an isolated sphere 
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on which small-scale fluctuations corresponding to a 
CDM power spectrum are superimposed. Here, we use 
Bertschinger's software COSMICS (Bertschinger 1995) in 
generating initial fluctuations. To incorporate the effects 
of fluctuations with longer wavelengths, the density of 
the sphere has been enhanced and a rigid rotation corre- 
sponding to a spin parameter, A, has been added. The 
initial condition of this model is determined by four pa- 
rameters A, Mtot, 08, in, and z c : the spin parameter is 
defined by 



GM. 



5/2 



(23) 



tot 



where J is the total angular momentum of the system, E 
is the binding energy when the system turns around, and 
M tot is the total mass of this sphere, which is composed 
of dark matter and gas; <7s.in is an rms mass fluctuation in 
a sphere of radius 8 h^ 1 Mpc, which normalizes the am- 
plitude of the CDM power spectrum; z c is the expected 
collapse redshift. If the top-hat density perturbation has 
an amplitude of Si at the initial redshift, Zi, we have 
z c = 0.36^(1 +Zi)-1 (e.g., Padmanabhan 1993). Thus, 
when z c is given, Si at Zi is determined. 

Our simulations assume a flat universe (fl = 1) with 
a baryon fraction of Slf, = 0.1 and a Hubble constant of 
H = 50 km s _1 Mpc . In this paper, we focus on the 
evolution of a top-hat density perturbation with a spin 
parameter A = 0.02. This is close to the minimum value 
of the spin parameter possible in a CDM universe, ac- 
cording to the results of the N-body simulation (Barnes, 
Efstathiou 1987; Warren et al. 1992). Other parameters 
are M tot = 8x 10 11 M©, cr 8 ,in = 0.5, and z c = 2. Then, 
the radius of the proto-galactic sphere is 32.17 kpc (in 
real space) at the initial redshift of Zi = 40. 

We carry out two simulations with the different num- 
bers of particles, N p . One simulation has N p = 9171 
(Model 1) and another has N p = 3071 (Model 2). In 
Model 1, the initial redshift is set to be Zi = 40, while 
in Model 2, Zi = 25. The numerical simulation should 
start from a redshift at which small scale perturbations 
are weak enough. The simulation with large N p can in- 
clude smaller-scale perturbations than that with small 
N p , and the perturbations on a smaller scale have a 
larger amplitude than large-scale perturbations in the 
CDM power spectrum. This is why the initial redshift 
of Model 1 is taken to be higher. In Model 1 (Model 
2), the masses of gas, star, and dark matter particles 
are 8.72 x 10 6 (2.61 x 10 7 ), 8.72 x 10 6 (2.61 x 10 7 ), and 
7.85 x 10 7 (2.34 x 10 8 ) M Q , respectively and the soft- 
ening lengths of gas, star, and dark matter particles are 
1.20 (1.73), 1.20 (1.73), and 2.50 (3.61) kpc, respectively. 



4. Results 

We simulate the evolution of the galaxy based on the 
two models from the redshift z% to the present, z = 0. 
The angular momentum of the whole system is conserved 
with sufficient accuracy: the error is less than 0.3 (1) % 
in Model 1 (Model 2). Figure 1 shows the evolution of 
the system in Model 1. These panels show the projection 
of the particles to the x-z plane, where we take the z- 
axis to be the initial rotational axis. At redshift of z = 
4.88, the system is expanding with the Hubble Flow. The 
amplitude of the initial top-hat density perturbation is 
chosen so that the system should turn around at z = 3.71. 
At z = 3.54, small clumps are formed and star formation 
begins in these clumps. These clumps have merged and 
the whole system has collapsed at z = 1.89. The stars 
are formed in a burst and gas particles diminish, changing 
into star particles. At z = 0.84, the system has already 
relaxed and is evolving almost passively. 

The history of the SFR is shown in figure 2. Figure 2 
also shows the SFR in Model 2. In spite of the difference 
in the number of particles, the two simulations yielded 
similar results. From figure 2, we see that the major 
star formation begins due to the collapse of whole sys- 
tem and continues for about 1 Gyr with a rate of about 
50 Mq yr _1 . At the age of universe t = 4 Gyr (z <~ 1.2), 
the star formation is almost quenched owing to the ex- 
haustion of the gas. 

The distribution of star particles at z = is shown 
in figure 3. Figure 3 shows the x-y and x-z projection. 
Both projections look quite similar. The final stellar sys- 
tem seems to be nearly spherical due to the low initial 
spin, in contrast to the results of Katz (1992) who made 
the disk system contracted along the rotational axis due 
to the high initial spin. 

Figure 4 shows the surface mass density profile of the 
final stellar system at z = in Model 1 and Model 2. 
These surface densities are obtained for the cylindrical 
shells of various radii in the x-y projection. We have 
confirmed that these profiles do not change in other pro- 
jections. The surface density profile in figure 4 is in good 
agreement with the de Vaucouleurs law, 



S = £ c cxp{-7.67 [(i?/i?ha„) 1/4 - l] } 



(24) 



where i? n ,m is the radius containing half the total mass 
and S e is the mass density at i? n ,m- These parameters 
are summarized in table 1. Here, -R n ,m is comparable to 
the softening length especially in Model 2. In order to 
examine the effect of the softening length of star particles, 
e s , which equals that of gas particles, we also carry out 
five simulations (N p — 3071) with the softening length 
e s = 0.25e s ,o, 0.5e Sj o, 0.75e s ,o, 1.25e S; o, and 1.5e S; o, where 
e s .o = 1.73 kpc denotes the standard value used in Model 
2. Figure 5 shows i? n ,m at z = as a function of e s . It 
is clear that e Si0 is an appropriate value, because i?h,m 
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Fig. 1. Time evolution of the system in Model 1. The upper, middle, and lower panels show the distributions of the dark matter, the 
gas, and the stars, respectively. Each panel measures 200 kpc across and shows the x-z projection of the particles, where we set the 
z-axis to be the initial rotational axis. 




Time (Gyr) Time (Gyr) 

Fig. 2. Star formation rate in Mq yr — 1 as a function of the age of universe in Model 1 (left) and Model 2 (right), respectively. 
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Y 



X 



X 



Fig. 3. The particle distribution of stars at z = 0. The left (right) panel shows the x-y (x-z) projection, where the z-axis is the initial 
rotational axis. Each panel measures 100 kpc across. 




Fig. 4. The surface mass density profiles in Model 1 (left) and Model 2 (right), respectively. The solid line shows the de Vaucoulcurs law 
given by equation (24) with parameters being listed in table 1. 
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Softening Length (kpc) 



Fig. 5. Half-mass radius as a function of the softening length of the 
star particles in the simulations with 7V P = 3071. A filled circle 
corresponds to the softening length in Model 2 (e B ,o = 1-73 kpc). 




Radius (kpc) 



does not depend on the softening length around £ Si o- It is 
worth noting that, in the range e s < 0.5e s ,o, the softening 
length is so small that the artificial two-body relaxation 
causes the rise of i?h,m- Figure 6 shows this effect clearly. 
In the case e s — e a .o, Rh,m does not change from z = 
0.58, at which the system has already relaxed, to the 
present epoch z = 0. In the case e s = 0.25e S! o, however, 
the central concentration is diluted and i?h,m goes up at 
z — 0. Thus, the resulting value of i?h,m appears to be 
trustful only over a limited range of e s , which includes 
e Si o adopted in Model 2. 

It is known that a violent relaxation for sufficiently 
deep central potentials reproduces the de Vaucouleurs 
law (Hjorth, Madsen 1991 and references therein). The 
deep potential causes the strong scattering which tightly 
bound particles experience near the center, and leads to 
the outer envelope of the de Vaucouleurs profile. In dis- 
sipationless collapse, a cold collapse is preferred to yield 
a deep enough central potential (e.g., van Albada 1982). 
However, the dissipation is a crucial component of galaxy 
formation. Carlberg et al. (1986) found that dissipation, 
leading to a deep central potential, could help get a good 



Fig. 6. The surface mass density profiles at z = 0.58 (short-dashed 
lines) and at z = (solid lines), respectively. The thick and normal 
lines correspond to the case of e s = e Sj o and the case of e s = 
0.25e Sl o, respectively. The arrows indicate the position of the half- 
mass radius. 
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de Vaucouleurs profile. We can infer that our numeri- 
cal simulations also result in the dc Vaucouleurs profile 
owing to the collapse combined with the dissipation. 

In Model 1, we analyze the kinematic properties, the 
rotation velocity and the velocity dispersion. The maxi- 
mum rotation velocity V m , is 25 km/s at the radius of 9 
kpc where the rotation curve flattens out, and the mean 
velocity dispersion, a, within the half-mass radius is 116 
km/s. V m /a is a important value characterizing the kine- 
matic properties of elliptical galaxies. V m /a = 0.2 means 
that this system is slowly rotating and supported by the 
velocity dispersion due to the low initial spin. This value 
resembles that of bright elliptical galaxies (e.g., Davies 
et al. 1983). 

Our simulation includes the chemical evolution pro- 
cess. As a result, each star particle has its own age 
and metallicity. From this information, we can obtain 
the SED (Spectral Energy Distribution) using the SSPs 
(Single Stellar Populations) in Kodama, Arimoto (1997). 
Consequently, we can analyze distributions of the lumi- 
nosity and color. Figure 7 shows the surface brightness 
(Lq.v pc~ 2 ) in V band as a function of the radius in 
Model 1 and Model 2. The surface brightness profile also 
matches the de Vaucouleurs law very well, using the half- 
light radius i?h,i and the luminosity I a at i?h,i given in 
table 1. 

Based on the luminosity of each star particle, wc also 
obtain the radial profile of the metallicity and color 
weighted by the luminosity. Figure 8 shows the pro- 
jected metallicity and the V — K color profile in Model 1. 
The metallicity has a gradient of A log Z/A log r ~ —0.3, 
which is consistent with the value A log Z/A log r — 
-0.2 ± 0.1 observed by Davies et al. (1993). The V - K 
color also has a similar amount of gradient which is 
caused mainly by the spatial change of metallicity. 

In table 1, we summarize the photometric properties 
of this stellar system at z = 0. The eighth column in 
table 1 gives the V — K color within a 10 kpc aperture, 
which is the same aperture as adopted in the observa- 
tion by Bower et al. (1992), who investigated the color- 
magnitude relation of early-type galaxies in the Coma 
and Virgo clusters. The half-light radius, the metallicity 
gradient, and the V — K color obtained in the simula- 
tion agree with those values observed for a My ~ — 20 
elliptical galaxy (e.g., Faber et al. 1989; Peletier et al. 
1990; Bender et al. 1992; Bower et al. 1992; Davies et al. 
1993). There is no significant difference between the two 
simulations which have different N p . Therefore, these 
simulations have a sufficiently large particle number and, 
hence, sufficient mass resolution to reproduce the main 
properties of the elliptical galaxy with enough accuracy. 



5. Summary and Discussion 

We simulated the galaxy formation from a top-hat den- 
sity perturbation with a small spin parameter, A = 0.02. 
The stellar system formed in the simulation is spherical 
and has a mass distribution which is in good agreement 
with the de Vaucouleurs law. This system was formed 
through turning around from the Hubble expansion and 
subsequent merging of small clumps originating in initial 
CDM density perturbations. We also obtained photo- 
metric properties by simulating star formation process 
which follows chemical evolution. These properties were 
found to be similar to the observed properties of ellipti- 
cal galaxies. The results of our simulation clearly show 
that, in a CDM universe, the proto-galaxy which has a 
spin-parameter as small as 0.02 evolves into an elliptical 
galaxy. 

Using three-dimensional SPH numerical simulation, 
Mori et al. (1999) investigated elliptical galaxy forma- 
tion. They assumed a virialized proto-galactic gaseous 
cloud as the initial condition. As Katz & Gunn (1991) 
showed, however, the gas never gets heated to the virial 
temperature in the hierarchical galaxy formation which 
was adopted in our simulations. The present simulations 
thus improve on their non-cosmological simulations. 

In our numerical simulation, the luminosity profile is 
in good agreement with the de Vaucloulcurs law as well 
as the mass distribution, in spite of the difference be- 
tween the half mass and light radii. Figure 9 shows the 
projected stellar mass-to-luminosity ratio in V band as 
a function of the radius in Model 1. Due to the metal- 
licity gradient, the mass-to-luminosity ratio decreases as 
the radius increases. This gradient of mass-to-luminosity 
ratio leads to the half-light radius larger than the half- 
mass radius. However, the de Vauclouleurs law is held 
for the luminosity, if the mass-to-luminosity ratio varies 
with respect to the radius as follows, 



E(r) 
J(r) 



= ^-exp 



-7.67 {R-^ - iC 4 ) - 1/4 ] • (25) 



In figure 9, the solid line corresponds to this equation 
with parameters of table 1. The simulation results match 
this line. Thus, this gradient of mass-to-luminosity ratio 
does not make the luminosity profile fail to reproduce the 
de Vaucouleurs law. 

It is notable that we obtain the projected metallicity 
gradient which is consistent with the observed one for 
elliptical galaxies. This gradient is caused by different 
star formation history at different sites. Figure 10 shows 
the time variations of SFR in the regions of inside and 
outside the half-light radius in Model 1. It is clear that 
the star formation holds on longer time in the inner re- 
gion than in the outer region. Since the residual gas is 
polluted by the past star formation, metal-rich stars are 
formed in the central region where the duration of the 
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Fig. 7. The surface brightness profiles in Model 1 (left) and Model 2 (right), respectively. The solid line shows the de Vaucoulcurs law 
with parameters being listed in table 1. 




Fig. 8. The projected metallicity (left) and the V — K color (right) versus the radius in Model 1. In the left panel, the solid line has the 
gradient of AlogZ/Alogr = —0.3. 
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Table 1. de Vaucouleurs law parameters and photometric properties. 





N P 




-Rh,m 




Rh.i 


M v 


V-K 






(Mq PC" 2 ) 


(kpc) 


(L QiV pc~ 2 ) 


(kpc) 




(< lOkpc) 


Model 1 


9171 
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4.49 


-20.34 


3.164 


Model 2 
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Fig. 10. The history of the star formation rate (SFR) in different 
regions in Model 1. The solid line is the SFR within the half-light 
radius, while the dotted line is that outside the half-light radius. 



Fig. 9. The projected stellar mass-to-luminosity ratio in V band as 
a function of the radius in Model 1. The solid line corresponds to 
equation (25) with parameters of table 1. 



star formation is long. 

Since the cooling of the gas overcame the feedback from 
massive stars, our simulations did not lead to outflow of 
the gas in the form of a galactic wind, which is often pos- 
tulated to explain photometric properties of the elliptical 
galaxy population (e.g., Larson 1974b, Arimoto, Yoshii 
1987). Although our numerical model ignored some im- 
portant feedback processes like the stellar wind and the 
Type la supernova, the system which we simulated might 
be too massive to allow a galactic wind. To confirm this 
conjecture, we need to simulate the systems of different 
masses. The ultimate goal of such an extensive attempt 
is to reproduce the color-magnitude relation of elliptical 
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galaxies (e.g., Bower et al. 1992). A forthcoming paper 
will discuss the properties of model elliptical galaxies as 
a function of galaxy mass. 
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